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Abstract 

We propose a new numerical approach to compute nonclassical 
solutions to hyperbolic conservation laws. The class of finite differ- 
ence schemes presented here is fully conservative and keep nonclassical 
shock waves as sharp interfaces, contrary to standard finite difference 
schemes. The main challenge is to achieve, at the discretization level, 
a consistency property with respect to a prescribed kinetic relation. 
The latter is required for the selection of physically meaningful non- 
classical shocks. Our method is based on a reconstruction technique 
performed in each computational cell that may contain a nonclassical 
shock. To validate this approach, we establish several consistency and 
stability properties, and we perform careful numerical experiments. 
The convergence of the algorithm toward the physically meaningful 
solutions selected by a kinetic relation is demonstrated numerically for 
several test cases, including concave-convex as well as convex-concave 
flux-functions. 
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1 Introduction 



State of the art 

We are interested here in the challenging issue of numerically computing non- 
classical solutions (containing undercompressive shocks) to nonlinear hyper- 
bolic conservation laws. Nonclassical solutions have the distinctive feature of 
being dynamically driven by small-scale effects such as diffusion, dispersion, 
and other high-order phenomena. Their selection requires an additional 
jump relation, called a kinetic relation, and introduced in the context of 
phase transition dynamics [281 E21 EOJ EU [H El [H] EQl H3j HH [261 EZ] , and 
investigated by LeFloch and collaborators in the context of general hyper- 
bolic systems of conservation laws (see [21] for a review) . 

From pioneering work by Hayes and LeFloch \13\ [H] it is now recog- 
nized that standard finite difference schemes do not converge to nonclassical 
solutions selected by the prescribed kinetic function. In fact, kinetic func- 
tions can be associated not only with continuous models, but with the finite 
difference schemes themselves. Achieving a good agreement between the 
continuous and the numerical kinetic functions has been found to be very 
challenging. 

In the present paper, we will show how to enforce the validity of the 
kinetic relation at the numerical level, and we design a fully conservative 
scheme which combines the advantages of standard finite differences and 
Glimm-type (see below) approaches. 

Nonclassical shocks and other phase transitions are naturally present 
in many models of continuum physics, especially in the modeling of real 
fluids governed by complex equations of state. This is the case, for in- 
stance, of models describing the dynamics of liquid- vapor phase transitions 
in compressible fluids, or of solid-solid phase transformations in materials 
such as memory alloys. For numerical work in this direction we refer to 

pg QE1 BUSIES]. 

Setting for this paper 

We restrict here attention to scalar conservation laws 

dtu + d x f{u) = 0, u(x,t)eR, (i,t)eRxI+, 
u(x,0) = u (x), 

and postpone the discussion of systems of conservation laws to the follow- 
up paper [3]. The above equation must be supplemented with an entropy 
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inequality of the form 

d t U{u) + 8 x F(u) < 0. (2) 

Here, t denotes the time variable, x the (one-dimensional) space variable, 
/ : R — ► M the flux function, and (U, F) is any strictly convex mathematical 
entropy pair. That is, U : M — > M is strictly convex and F : R — > R is given 
by 7 1 ' = 17'/'. Equations (fTl) and Q are imposed in the distributional sense. 

We rely here on the theory of nonclassical solutions based on kinetic re- 
lations, established in [21 . The flux / is assumed to be nonconvex, which is 
the source of mathematical and numerical difficulties. From the mathemat- 
ical standpoint, a single entropy inequality like ([2J does not suffice to select 
a unique solution. This can be seen already at the level of the Riemann 
problem, corresponding to ([TJ-Q when uq has the piecewise constant form 

/ \ f ui, x < 0, , . 

U0(X) = {u r , x > 0, (3) 

ui and u r being constant states. The Riemann problem admits (up to) a 
one-parameter family of solutions (see Chapter 2 in [21 J. However, these 
solutions contain discontinuities violating the standard Lax shock inequali- 
ties, which are referred to as nonclassical. They are essential from the phys- 
ical standpoint, and should be retained. This non-uniqueness can be fixed 
however, provided an additional algebraic condition, the so-called kinetic 
relation, is imposed on each nonclassical shock. Consider a shock connect- 
ing a left-hand state U- to a right-hand state n+ and propagating with the 
speed o given by the usual Rankine-Hugoniot relation, that is, 

, s f u-, x<at, , > f(u + )-f(uJ) 

u(x,tj = < , a = a(u-,u + ) = . (4) 

v ' \ u+, x > at, y + ' u+-u- K ' 

The kinetic relation takes the form 

n+ = ip (it-) for all nonclassical shocks, (5) 

where ^ is the so-called kinetic function. Equivalently, denoting by ip ^ the 
inverse of the kinetic function it may be preferable to write n_ = f~'"(u + ). 
The kinetic relation implies that the right-hand (respectively left-hand) state 
is no longer free (as in a classical shock wave) but depends explicitly on the 
left-hand (respectively right-hand) state. 
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Objectives in this paper 

At the numerical level, several strategies exist in the literature in order to 
take into account the kinetic relation ([5j . We can distinguish between diffuse 
interface methods and sharp interface methods. 

In the first approach, one assumes that the kinetic relation is derived 
from an augmented continuous model and, in order to take into account the 
internal structure of nonclassical discontinuities, one attempts to resolve the 
effects dues to (small) diffusive and dispersive terms that generate them. It 
is then possible to construct conservative schemes that mimic at the numer- 
ical level the effect of the regularized models. Due to the great sensitivity of 
nonclassical solutions with respect to small scales and numerical diffusion, 
it turns out that numerical results are satisfactory for shocks with moderate 
amplitude, but discrepancies between the exact and the numerical kinetic 
function arise with shocks with large amplitudes and in long-time computa- 
tions. For this circle of ideas we refer the reader to |13l fT4"] , and the follow-up 
papers |22HZ||g|. 

In the second approach, small scale features are not explicitly taken into 
account. Instead, the kinetic relation is included, in a way or another, in 
the design of the numerical scheme. This is the case of the random choice 
and front tracking schemes. It should be mentioned here that the Glimm 
scheme and front tracking schemes do converge to exact solutions even in 
presence of nonconclassical shocks; see \20\ |2T| [23] for the theoretical aspects 
and Chalons and LeFloch [9 for a numerical study of the Glimm scheme. 
These schemes require the explicit knowledge of the underlying nonclassical 
Riemann solver, which may be expensive numerically, and this motivated 
the introduction of the so-called transport-equilibrium scheme by Chalons 

In [16] , Hou, LeFloch, and Zhong proposed a class of converging schemes 
for the computation of propagating solid-solid phase boundaries. More re- 
cently, Merckle and Rohde developed a ghost-fluid type algorithm for a 
model of dynamics of phase transition. These schemes provide satisfactory 
numerical results, as nonclassical discontinuities are sharply and accurately 
computed. Although the convergence of the methods was demonstrated 
numerically, their main drawback in practice is similar to the Glimm-type 
schemes and the property of strict conservation of the conservative variable 
u fails. 

Building on these previous works, our objective in this paper is to de- 
sign a fully conservative, finite difference scheme for the approximation of 
nonclassical solutions to the hyperbolic conservation law ([T]). Our basic 



4 



strategy relies on the discontinuous reconstruction technique proposed re- 
cently in Lagoutiere |18| [T§] which has been found to be particularly efficient 
to computing classical solutions of ([I]) with moderate numerical diffusion. 

In our approach below, the kinetic function ip b is included explicitly in 
the algorithm, in such a way that nonclassical shocks are computed (essen- 
tially) exactly while classical shocks suffer moderate numerical diffusion. To 
validate our strategy we perform various numerical experiments and, in par- 
ticular, draw the kinetic function associated with our scheme. As the mesh is 
refined, we observe that the approximate kinetic function converges toward 
the analytic kinetic function. The scheme also enjoys several fundamental 
stability properties of consistency with the conservative form of the equation 
and (like the Glimm scheme) with single nonclassical discontinuities. 



2 Nonclassical Riemann solver with kinetics 
Assumption on the flux-function 

We describe here the nonclassical Riemann solver introduced and investi- 
gated in LeFloch [21 . Note in passing that this solver was later extended 
in [23] to include also a nucleation criterion. 

Consider the problem ([T|-([2|-([5]) for a given Riemann initial data ([3]). 
Throughout this paper we assume that the flux / is either concave-convex 
or convex-concave, that is, satisfies the conditions (for all 

n/»>0, f'(0)/0, lim M _ +00 /'(«) = +<», (6) 

or 

u/"(u)<0, f"(0)/0, lim M _ +00 /'(«) = -oo, (7) 

respectively. The functions f(u) = v? + u and f(u) = —u 3 — u are proto- 
types of particular interest, used later in this paper for the validation of the 
proposed numerical strategy. 

Let ip* : R — ► R be the unique function defined by (p\0) = and for all 
u ^ 0, tp\u) ^ u is such that the line passing through the points (u, f(u)) 
and (tp\u), f(if^(u))) is tangent to the graph of / at point (ip^u), f(tp*(u))): 

k f{u)-f{y\u)) 

f kv w) = wr~\ — • 

u — <p\u) 

This function is smooth, monotone decreasing and onto thanks to (pj or 0. 
We denote by tp~^ : R — ► R its inverse function. 
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Concave-convex flux functions 

Let us assume that / obeys ^ and let 99^ : K — > M be a kinetic function, that 
is (by definition) a monotone decreasing and Lipschitz continuous mapping 
such that 

<Po(u) < ^( u ) < ^(u), u > 0, 

ip\u) < ip\u) < ifliu), u<0. U 

From y^, we define the function p^ : R — > WL such that the line passing 
through the points (u,f(u)) and (<p^(u), f(ip*(u))) with also cuts the 

graph of the flux function / at point (</?"(«), /(</?"(«))) with 7^ it and 

it — (/? b (n) it — <p$(u) 

The nonclassical Riemann solver associated with ([IJ-Q-Q-Q is given as 

follows. 

When ui > 0: 

(1) If u r > m, the solution is a rarefaction wave connecting ui to u r . 

(2) If ti r G [<p$(ui),Ui), the solution is a classical shock wave connecting 
ui to u r . 

(3) If it r S (<^ 1, ('U/), 9?" the solution contains a nonclassical shock 
connecting ui to cp (ui), followed by a classical shock connecting </3^(?xz) to 
u r . 

(4) If u r < (p (ui), the solution contains a nonclassical shock connecting 
to (p°(ui), followed by a rarefaction connecting tp (u{) to u r . 

When ui < 0: 

(1) If u r < u/, the solution is a rarefaction wave connecting u\ to it r . 

(2) If u r £ [uz, (?//)), the solution is a classical shock wave connecting 
ui to u r . 

(3) If it r S (tp$ (u[) , (p b (ui)) , the solution contains a nonclassical shock 
connecting ui to ip (ui), followed by a classical shock connecting <p (ui) to 
it r . 

(4) If u r > ^ the solution contains a nonclassical shock connecting 
ui to (p°(ui), followed by a rarefaction connecting (p (ui) to u r . 

Convex-concave flux functions 

We next assume that / satisfies the condition Q. Let <p b : R — > R be a 
kinetic function, that is, a monotone decreasing and Lipschitz continuous 
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map such that 

p b (u) < p\u) < <p-*(u), u<0, 

ip-\u) < p\u) < (fl(u), u>0. { ' 

We then define p(u, v ) G R if v ^ u and v ^ ip\u) by 

f(p(u,v))-f(u) = f(v)-f(u) 
p(u, v) — u v — u 

with p{u, v) ^ u and i>) ^ v, and extend the function p by continuity 
otherwise. Note that p$(u) = p(u, (f (u)) where cp$ is defined as in the case of 
a concave-convex flux function. The nonclassical Riemann solver associated 
with @-@-((3])-([5]) is given as follows. 
When ui > 0: 

(1) If u r > ui, the solution is a classical shock connecting u\ to u r . 

(2) If u r G [0,itj), the solution is a rarefaction wave connecting m to u r . 

(3) If u r G (y> (mz), 0), the solution contains a rarefaction wave connecting 
Ui to p~ b (u r ), followed by a nonclassical shock connecting p~ b {u r ) to ii r . 

(4) If u r < (itj), the solution contains: 

(i) a classical shock connecting uz to p (u r ), followed by a nonclas- 
sical shock connecting tp (u r ) to u r , if uz > p{<p~ (u r ), u r ). 

(ii) a classical shock connecting ^ to u r , if uz < p(p~ ] "(u r ), u r ). 
When m < 0: 

(1) If « r < ui, the solution is a classical shock connecting u\ to u r . 

(2) If u r G (ui, 0], the solution is a rarefaction wave connecting ti/ to u r . 

(3) If ii r G (0, (/^(uz)), the solution contains a rarefaction wave connecting 
Ui to ^"'(u,.), followed by a non classical shock connecting <p~ b (u r ) to u r . 

(4) If u r > (p^(ui), the solution contains: 

(i) a classical shock connecting uz to (/? (w r ), followed by a nonclas- 
sical shock connecting p (u r ) to tz r , if ui < p(<p (u r ) , u r ) . 

(ii) a classical shock connecting ui to u r , if uz > p(p~ ] "(u r ), u r ). 
Observe that the convex-concave case can in principle be deduced from 

the concave-convex case, by replacing / by — / and x by —x. Nevertheless, 
it is useful to keep the above two descriptions in mind, since there is a 
dramatic difference between the Riemann solvers: the nonclassical shock 
always connects ui to p?(ui) in the concave-convex case, and p~ b (u r ) to u r 
in the convex-concave case. The numerical method we are going to describe 
must take this feature into account, and as we will explain it is necessary to 
take into account both <p b and in the design of the scheme. 
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3 Motivations and difficulties 



Notation 

Our aim is to design a scheme for the numerical approximation of the non- 
classical solutions to ([T])-(|2])-(|5|. To this end, we consider the general class of 
finite volume methods. Introducing constant space and time lengths Ax and 
At for the space and time discretization, we can set Xj+1/2 = jAx, j G Z, 
and t n = nAt, n E N. The discretization consists, at each time t n , of a piece- 
wise constant function x i— ► u u (x, t n ) which should be an approximation of 
the exact solution u(x, t n ) on the cell Cj = [xj_i/ 2 ; x J+1 / 2 ): 

u u (x,t n )=Uj, xGCj, jeZ, n £ N. 

Here, v refers to the ratio At/Ax. The initial data at the time t = is 
denoted by uq and we define the sequence (u^)j & z- 

1 f x j+l/2 



= — / u (x)dx, j G Z. (10) 

X ^ x j-l/2 

The starting point in the conception of our algorithm is a few conven- 
tional interpretation of the constant values uj, j £ Z. As suggested by 



the proposed initialization (10), u™ is usually, and rightly, seen as an ap- 
proximate value of the average on cell Cj of the exact solution at time t n . 
Integrating equation ([!]) over the slab Cj x [f n ,t n+1 ] and using Green's for- 
mula, it is thus natural to define {u 7 j +1 )j from (iij)j and a conservative 
scheme of the following form 

< +1 = «? - - tf-i/a)' ( u ) 

where /JYi/2 re P resen ts an approximate value of the flux that passes through 
the interface Xj + \/2 between the times t n and t n+1 . 

Here, we shall also consider vJ- as a given information, on cell Cj and at 
time t n , on the structure of the exact Riemann solution associated with 
inital states u\ = u 1 j_ 1 and u r = which will develop at the next times 
t > t n . At this stage, one easily realize that if this information is precise 
(i.e. close to what will really happen), then we should be be in a good 
position to define accurately the numerical fluxes and then predict 



the approximate values of the solution at time t 



n+l 
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Linear advection equation 

As a first illustration, let us consider the linear advection with constant ve- 
locity a > 0, that is, the scalar conservation law with flux f(u) = au. In 
this case, the weak solution to the initial-value problem for ([!]) is unique, 
and is given explicitly as u(t,x) = uq(x — at). Hence, neither the entropy 
condition ^ nor the kinetic condition ^ are necessary. The basic scheme 
for approximating this solution is the so-called upwind scheme and corre- 
sponds to the choice /jYi/2 = au j ^ or an 3 ^ R ecan that the CFL 
condition aAt/Ax < a for a given a < 1 is mandatory for the stability 
of the procedure. Figure [T] (left-hand) shows the corresponding numerical 
solution at time t = 0.25 for a = 1, a = 0.5 and u\ = 1, u r = in ([3]). The 
mesh contains 100 points per unit. We observe that the numerical solution 
presents a good agreement with the exact one but contains numerical diffu- 
sion. We propose the following interpretation. In some sense, the value iff 
that we consider as an information on the Riemann solution associated with 
initial states u\ = u 1 j_ 1 and u r = Uj +1 is sufficient to correctly approach 
this solution when defining /JY1/2 = aM i > no ^ enou gh to avoid the nu- 
merical diffusion. Note that the latter is expected but not hoped. In the 
present situation, the fact is that we actually know what will happen in the 
future, namely a propagation of the Riemann initial states (ui = and 
u r = Uj + i) with speed a. In particular, no value different from u 1 j_ 1 and 
is created so that information given by it" is clearly not optimal. In the 
process of calculation of the numerical flux /™ + i/2' we are ^ us tempted to 
add more information in the cell Cj when replacing, as soon as possible, the 
constant state it" with a discontinuity separating on the left and Uj +1 
on the right, and located at point Xj G Cj. In the forthcoming developments, 
the left and right states of this reconstructed discontinuity will be noted 
and Uj r , respectively. Hence, we have here 

= u% r = u? +1 . (12) 

We claim that this provides better information for calculating /j^]/ 2 than 
the original one. Such a reconstruction is due to conserve u in order to be 
relevant, which defines Wj by the following constraint 



(Xj - Xj-.x^Ujj + (x j+1 / 2 - Xj)u] >r = (x j+1/2 - Xj-xj^U 

which equivalently recast as 



^~^Ax. (13) 

, n _ „.n \ I 

U j,r U j,l 
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Figure 1: Linear advection - upwind scheme (left-hand) and reconstruction 
scheme (right-hand). 
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Then, the reconstruction is possible provided we have < < 1, with 



j 

1i A <r- W A 

a 3 ~ ,.n _ ,.n • V i4 J 



Now, let us introduce A^ +1 / 2 the time needed by the reconstructed discon- 
tinuity to reach the interface £7+1/2 (recall that a > 0). We clearly have 



1 - d n 
At j+1/2 = L Ax. 

In this case, the flux that passes through £7+1/2 between times t n and t n+1 = 
t n + At equals f(u^ r ) until t n + At j+1/2 , and after (if At j+1 / 2 < At). 

Therefore, we propose to set now 



Ai/jVi /2 = min(At i+1/2 , At)/(^ r ) + max(At - At J+1/2 , 0)/(^ ; 



On Figure [T] (right-hand) , we have plotted the numerical solution given by 
this new numerical flux, leading to the so-called reconstruction scheme. The 
parameters of the simulation are the same than those of Figure [l] (left-hand). 
We see that the more precise informations we have brought on each cell Cj 
for calculating the numerical fluxes make the scheme less diffusive than the 
original one. This strategy was proposed (and is discussed in further details) 
in |18| [T§] (see also |1U| [T7]). In particular, it is shown therein that the 
numerical solution presented in Figure [T] (right-hand) is exact in the sense 
that equals the average of the exact solution on Cj , that is 



1 r x i+i/2 

u] = — / u(x,t n )dx, j € Z, n G N. (15) 



The corresponding numerical discontinuity separating ui and u r in then 
diffused on one cell at most. 



Godunov scheme with a nonclassical Riemann solver 

As a second illustration, let us go back to the problem (|l])-(|2|-([5| with a 
general concave-convex (or convex-concave) flux function / with however, 
for the sake of clarity, 

f'(u) > 0, ueR. (16) 

Here, we focus ourselves on a particular Riemann initial data ^ such that 
u r = (p (ui). In other words, the kinetic criterion is imposed on the initial 
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discontinuity. The exact solution then corresponds to the propagation of this 
discontinuity with speed a(ui,u r ) > given by Rankine-Hugoniot relation: 



a(u h u r 



f(Ur) ~ f(ui) 
U r - Ul 



(17) 



Figure [2] (left-hand) represents the numerical solution given by the upwind 
scheme /" + i/2 = ) a * time t = 0.1, for f(u) = + u and ui = 1. The 
kinetic function is taken to be <p^(u) = — 0.75 u so that u r = —0.75. 




Figure 2: Propagating nonclassical shock - upwind scheme (left-hand) and 
reconstruction scheme (right-hand). 



We observe a strong disagreement between the numerical solution and 
the exact one. Indeed, the former is made of a (classical) shock followed 
by a rarefaction wave while the latter is a single (nonclassical) shock from 
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ui to u r . It is then clear that the usual upwind scheme (as many others 
actually) is not adapted for the computation of nonclassical solutions. The 
next result states that the upwind scheme always converges towards the 
classical solution of ([TJ-Q. This scheme is then adapted for the computation 
of classical solutions only. 

Property. Assume that uq £ L°°(M) and f is a smooth function satis- 



fying (16). Then, under the CFL condition 

^max|/»|<l, 
Ax ueA 

with A := [mm x uo(x) , m&x x uq(x)] the upwind conservative scheme 
with /" + i/2 = f( u ^) converges towards the unique classical solution of 

To establish this property, we only need to observe that, under the as- 



sumption (16) (propagation is only in one direction), the upwind scheme 
is equivalent to the standard Godunov scheme associated with the classical 
Riemann solver of ([T])-(|2j Then, standard compactness and consistency ar- 
guments apply and allow us to conclude that the scheme converges towards 
the unique classical solution. 

Obviously, the above property also holds if / is assumed to be decreasing 
if we define /™ + i/ 2 = f( u j+i)- 

4 A conservative scheme for nonclassical entropy 
solutions 

Preliminaries 

In view of the discussion in the previous section and in order to better eval- 
uate the numerical fluxes /" +1 / 2 , let us obtain some information beyond u™ 
on cell Cj. In the present instance of an isolated propagating discontinuity, 
it is expected that the Riemann solution associated with initial states t*"_i 
and Uj + i simply propagates the initial discontinuity. This is actually true if 
itj_i = v-t and Uj +1 = (p b (ui), or more generally if u^ +1 = yf(u'j_ 1 ). So that 
here again, we propose to replace the constant state with a discontinuity 
separating and u^ r and located at point Xj given by (13), as soon as 
possible i.e. when < d™ < 1. We take 

u^ = ^« +1 ) and ul r = v\u™_ x ). (18) 
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Note that this reconstruction is equivalent to (12) provided that u 1 J_ 1 = u\ 
and u!L_i = (f b (ui), or more generally u r - +1 =ipP(u r j"_ 1 ). Then, under the 



assumption (16), we again naturally set 



Atf, 



i+1/2 



min(At j+1/2 , At)f(u% r ) + max(Ai - At i+1/2 , 0)/(u£,) 



with now 



At 



i+1/2 



<t{u" 



-Ax. 



(19) 



Figure [2] (right-hand) highlights the benefit of such a reconstruction. The 
numerical solution now fully agrees with the exact one and is moreover free 
of numerical diffusion (the profile is composed of a single point). We will 



show below that it is exact in this case, in the sense that (15) is still valid 
as in the linear case. 



The scheme 

On the basis of the above motivations and illustrations, we follow the de- 
scription of our algorithm by considering the general situation. Assuming 
as given a sequence (uj)j^z at time t n , it is thus a question of defining its 
evolution towards the next time level t n+ . More precisely, and in the con- 
text of a finite volume conservative scheme, we have to define the numerical 
fluxes (/™ +1 / 2 )jez coming in ( JTTj ) . For that, we still assume 

either f(u) > for all u, or f'(u) < for all u, (20) 

so that propagation is in one direction only. According to the previous 
section, information in cell Cj is understood as an element of the inner 
structure of the Riemann problem associated with initial states u^Lj and 
u j"+l- This one will be used to compute either 7™ + i/2 f'( u ) — 0) or fj-x/2 

uf ./";»)'<>). 

In Section [2] it is stated that the Riemann problem associated with initial 
states Uj_i and Uj +1 may contain a nonclassical shock between and 
cp if the function is concave-convex (and between ip~^{u^ +1 ) and u™ +l 

if the function is convex-concave) . 

Recall that these nonclassical waves are difficult to capture numerically 
and require special attention. (We have shown in the previous section that as 
many others, the upwind scheme does not suit.) Instead of considering it" as 
a sufficiently accurate information for the structure of the Riemann solution 
associated with the initial states u"_ 1 and it^ +1 , we propose to replace it 
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(whenever possible) with a discontinuity separating = <p t, (u™ +1 ) on the 
left and u™ r = ip (it™^) on the right, and located at point Xj 6 Cj. In other 
words, we propose to introduce in the cell Cj the right (respectively left) state 
cp^(u'j_ 1 ) (respectively ^^(ttj+i)) of the nonclassical discontinuity which is 
expected to be present in the Riemann solution associated with u'j_ 1 and 
Uj +1 (depending on if / obeys (|6j) or (j7|). As in the previous section, one 
requires the reconstructed discontinuity to satisfy the conservation property 



(13) and to be located inside Cj, that is < cf™ < 1 with given in (14). 



Here, we let u™ z 

Then, we naturally set for all j £ Z: 
(i) if / is non-decreasing 



iff if dj given in (14) does not belong to [0, 1]. 



At/? 



J+l/2 



' min(At i+1/2 ,At)/(up+ max(At - At j+1/ 
At/(^), 



2,0)/« 



< d% < 1, 



otherwise, 



with 



At 



i+1/2 



a ( u j,n u j,r) 



Ax. 



(21) 
(22) 



(zi) if / is non-increasing: 
min(At 



At/JL 



1/2 



i _ 1/2) At)/(^)+ 



max(At-At i _ 1/2 ,0)/(< 



j,rJi 



At/(«7), 



< d? < 1, 



otherwise, 



with 



At 



J- 1/2 



<7" 



-Ax. 



(23) 

< u iv u irr~' (24) 

Note that contrary to the linear advection (see the first illustration in the 
previous section), the local time step At J+1/ / 2 (respectively Atj_i/ 2 ) given 
by (22) (respectively (24)) is now only a prediction of the time needed by 
the reconstructed discontinuity to reach the interface 2^+1/2 (respectively 
Xj-i/2)- The prediction step is however exact in the case of an isolated non- 
classical discontinuity (see the second illustration in the previous section) 
and more generally as soon as and u r - +1 verify = (p"(u'j_ 1 ). 

Observe that the proposed scheme belongs to the class of five-point 



schemes, since depends on u 7 j_ 2 , 



-1j 



'j'+i 



and it? 



i+2- 
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Stability and consistency properties 



We now state and prove important properties enjoyed by our algorithm. 

We assume that the flux / satisfies the monotonicity condition ( |20[ ) and 
either the concave-convex or concave-convex conditions Q or ([7]) respec- 
tively. Then, under the CFL restriction 



At 
Ax 



max|/'(u)| < 1, 



(25) 



where the maximum is taken over all the u under consideration, the conser- 



vative scheme (11) with /jYi/2 defined for all j € Z by (21)-(23) is consistent 
with @-@-((5fm the following sense. 



Property 1 (Flux consistency.) Assume that u := u 
then f» +1/2 = f(u) iff>0 and f»_ 1/2 = f{u) if f < 0. 



Property 2 (Classical solutions.) Assume thatu r -_ 2 , Uj, an d 

-u™ +2 belong to the same region of convexity of f. Then the definition u™ +1 
given by the conservative scheme (11)-(21)-(2S) coincides with the one given 
by the usual upwind conservative scheme. Then it obeys all the usual stability 
properties provided by this scheme. In particular, the strategy is convergent 
if the whole discrete solution belongs to the same region of convexity of f. 



Property 3 (Isolated nonclassical shock waves.) Let u\ and u r be two 

initial states such that u r = ip b (ui). Assume that u° = U\ if j < and 
ifi = u r if j > 1. Then the conservative scheme (ll)-(21)-(23) provides an 



exact numerical solution on each cell Cj in the sense thai 



U ; 



1 

Ax 



u(x, t n )dx, 



j € Z, n e N, 



(26) 



-1/2 



where u denotes the exact Riemann solution of {^p-{§p-(^p-{5p given by 
u(x,t) = ui if x < o~(ui,u r )t and u(x,t) = u r otherwise, and is conver- 
gent towards u. In particular, the numerical discontinuity is diffused on one 
cell at most. 



The following comments are in order. Property (i) shows that the pro- 
posed numerical flux function is consistent in the classical sense of finite 
volume methods. Properties (ii) and (Hi) provide us with crucial stabil- 
ity/accuracy properties. They show that the method is actually convergent 
if the solution remains in the same region of convexity of / (see (ii)) or, 
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more importantly, the solution consists in an isolated nonclassical discon- 
tinuity satisfying the prescribed kinetic relation (see (Hi)). We emphasize 
that all of the conservative schemes proposed so far in the literature violate 
the latter property. 



Proof of Property [TJ (i) If u := 



u 



d", 



<P (u) 



then 



^(u) - (p-^u)' 



The property < d™ < 1 means min(<^ ^(u), ip"(u)) <u< max(y i P^( u )) 
and cannot hold, since u and ^(u) do not have the same sign for all u. Then, 
we obtain ff +1/2 = f(u) if /' > and /«_ 1/2 = f(u) if /' < by Q)-(|23). 

Proof of Property [2} Assume without restriction that /' > and recall 
that < dj_i < 1 and < d™ < 1 respectively means that 



min(^ b (n"),^«_ 2 )) < < max(^(«"),^K_ 2 )) 



and 



min(^(^ +1 )VK-i)) < u^max^^),^^)). 



These inequalities are not valid since by definition u and y (it) do not belong 
to the same region of convexity of /. By (21)-(23), the numerical fluxes 
fj±i/2 c °i nc ides with the usual upwind fluxes and the conclusion follows. 

Proof of Property [3[ First, note that there is no relevant reconstruction 
in the first iteration. Indeed, the property < < 1 reads as follows if 
j < or j > 1, 



< e£? < 1 if and only if 



min(y? b (ui),p b (ui)) < u t < max(y \ui), cp°(ui)), j < 0, 
min((p~''(u r ),if'(u r )) <u r < m&x.((p~ ] " (u r ) , (p"(u r )), j > 1, 



which again cannot hold (see (i) below), while if j = or j = 1, the relation 
u r = (f > (ui) and definition (14) give 



d'! 



U r 



- Uj 

■ ui 
u r 

■ ui 



1. 

0. 



0, 
1, 
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so that the reconstructions exist but are trivial: ui = <p~"(u r ) (respectively 
u r = (p b (ui)) takes the whole cell associated with j = (respectively j = 1). 
Assume now without restriction that / is non-decreasing and let At be 



such that (25) holds. After one time step At, the exact solution given by 
u(x, At) = ui if x < a(ui,u r )At and u(x, At) = u r otherwise is such that 



1 

Ax 



u(x, At)dx 



x i-l/2 



cr(ui,u r )^(u r 



ui), 



But recall that a(ui,u r ) is given by (17) so that we have 



Ax 



that is 



'j+l/2 



u(x, At)dx 



-1/2 



1 

Ax 



ui 
u r 

U r 

"j + l/2 



/(«i))> 

f(ui)), 

fM), 



3<0, 
3 = h 
3 > 1- 



3 <0, 
3 = h 
3 > h 



(27) 



u(x, At)dx, 



3 G 



(28) 



(29) 



-1/2 



The identity (26) is then proved for the first iterate. 



What happens now in the next time iteration ? At this stage, it is first clear 
(see the previous discussion just below) that only cell C\ is going to be dealt 
with a reconstruction. Now, the main point of the proof lies in the fact 
that the reconstructed discontinuity in this cell actually joins the expected 
states <£> (uq) = <p~^(u r ) = ui and ^(uq) = (f b (ui) = u r and is located 
exactly at point x = a(ui,u r )At by the conservation property (29). In other 
words, we have reconstructed the exact solution at time t = At. To derive 



the required identity (26) for the second iterate, it is sufficient to recall that 
by Green's formula the conservative scheme (11) with /" +1 / 2 defined for all 
j G Z by (21)-(23) is equivalent for n = 2 to average the evolution of this 
exact solution up to time t 2 = 2 At. And the process is going on in a similar 
way for the next time iterations, which proves the result. 



5 Numerical experiments 

We mostly consider here the flux f(u) := u 3 + u, thus / is concave-convex 
in the sense given in the second section. For the entropy-entropy flux pair 
(U, F) required in ([2]), we use 

U(u):=u 2 , F(u) :=^u 4 + u 2 . 
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Easy calculations lead to explicit formulas for tp^ and <p> 

<ft = <P~^ = -2it, tpo{u) = -u. 

Moreover, we have here ip^(u) = —u — p^{u). 

The choice of the kinetic function ipr must be in agreement with relations 
dsl) with (p^ and c/?q just calculated. Here, we will choose the kinetic function 

<p\u) = -(3u, 0€ [0.5,1), 

which, as observed in Bedjaoui and LeFloch [3], can be realized by an aug- 
mented model based on nonlinear diffusion and dispersion terms. In the 
following, we will take j3 = 0.75. 



Test A. Let us check Property[3]numerically, which is concerned with the 
exact capture of isolated nonclassical shocks. Thus, consider the following 
nonclassical shock as a Riemann initial condition 



u (x) 



4, x < 0, 
^(4) = -3, x > 0, 



The numerical solution shown in Figure [3] is exact everywhere but in the 
single cell containing the nonclassical shock. (We sometimes use a piecewise 
constant representation in the figure, in order for the interpretation of the 
numerical solutions to be easier.) However, as expected, the value in this cell 



coincides with the average of the corresponding exact solution (see ( 26 ) ) , and 
allows (after reconstruction) to recover the exact location of the discontinuity 
(using the conservation property of scheme) . This property explains why the 
numerical solution stays sharp when the time evolves. 



Test B. In our second test we consider the Riemann problem with initial 
data 

4, x < 0, 
-5, x > 0, 



u (x) 



whose solution is a nonclassical shock followed by a rarefaction wave. The 
two left-hand curves in Figure [4] are performed with Ax = 0.01 and Ax = 
0.002, respectively. The nonclassical shock, as previously, is localized in a 
single computational cell. 

The right-hand figure represents the logarithm of the L 1 -error (between the 
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exaci 
exact 
exact 

lumeri 
lumori 
lunvii 



initial data t=0.00 - 
1 solution at t=0.Q1 - 
:t solution at t=0.02 
:t solution at t=0.03 
I solution at t=0. 01 - 
I solution at t=0.02 - 
I solution at t=0.03 



Figure 3: Test A - Nonclassical shock - 30 points 



exact and the numerical solution) versus the logarithm of Ax. The numerical 
order of convergence is about 0.8374. 

Test C (Figure [5J. Now, we choose another Riemann initial condition 
which develops a nonclassical shock followed by a classical shock: 

( v / 4, x < 0, 
U0{X) = {-2, x > 0. 

We can make the same observation as previously, concerning the nonclassical 
shock; it is sharply captured and arises in a small spatial domain. However, 
note here that the classical shock does contain some numerical diffusion: in 
fact, our scheme is exactly the upwinding scheme if the values of the solution 
remains in a given convexity region for the flux /. 

Once again, the plot with the L 1 -error shows the numerical convergence 
with order about 0.9999. 

Test D (Figure [6]). We now take an initial data composed of two non- 
classical shocks that interact: 

r 4 = ip~\-3), x < 0.1 
Uo (x) = < -3, 0.1 < x < 0.2 

{ 2.25 = v? b (-3), x > 0.2. 

The computation is performed with Ax = 0.05 and plotted at four successive 
times t = 0, 0.0010, 0.0017, and 0.0020. We observe that the two nonclassical 
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att=0 - 

exact solution at t=0.03 - 
numerical solution - 100 points 
numerical solution - 500 points 




Figure 5: Test C - Nonclassical and classical shocks 
(log(E L i) versus log(Ax)) 



convergence 
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shocks cancel each other at the interaction, and generate a single classical 
shock, in accordance with the general theory in [21 . 



att=0.000 - 
exact solution at t=0. 010 - 
exact solution at t=0. 017 
exact solution at t=0.020 
numerical at 1=0.010 - 
numerical at 1=0.017 - 
numerical at t=0.020 



Figure 6: Test D - Interaction of two nonclassical shocks 



Test E (Figure Next, we consider the periodic initial condition 

u (x) = sin (J^J , 

with periodic boundary conditions u(— 0.5, t) = u(0.5, t). The exact solution 
is not known explicitly, so we compare our numerical solution with the 
solution generated by Glimm's random choice scheme [T2j in which we have 
replaced the classical solver by the nonclassical solver described in Section [2] 
We use here van der Corput's random sequence (a„), defined by 

m 
k=0 

where n = Y^k=Q^h^ \ ih £ {Oil}) denotes the binary expansion of the 
integer n. Figure [7] represents the solutions at the times t = 0, 0.25 and 
0.5 for our scheme with Ax = 0.01 and with Ax = 0.0001, and for the 
Glimm scheme with Ax = 0.0001 (to serve as a reference). The two meth- 
ods strongly agree. Roughly speaking, the increasing parts of uq evolve as 
rarefactions, while the decreasing parts are compressed and develop in a 
classical shock and, then, when left- and right-hand states at the shocks 
change sign, nonclassical shocks (which do satisfy the expected kinetic rela- 
tion) and new faster classical shocks on the right-hand side arise. 



23 



t=0 






-0.4 


-0.2 


0.2 


0.4 












t=0 










100 pts 


t=0.50 


1 








-000 pts 


t=0.50 " 










Glimrir. 


t=0.50 














0.5 



























\ 

\ 










0.5 
























-1 














-0.4 


-0.2 


0.2 


0.4 



Figure 7: Test E - Periodic initial data - reconstruction scheme and Glimm 
scheme 
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Test F (Figure [8J). To illustrate the behavior of convex- concave flux 
functions, we finally compute two Riemann solutions with opposite flux 
f(u) = —u 3 — u (so /' < and the solutions move from right to left) and 
the same kinetic function t£> (u) = —0.75 u: the first one (left-hand figure) 
corresponds to the initial data 

f -4, x<0, 
U0(X) = \ 4, *>0, 

and develops a rarefaction wave and a nonclassical shock; the second one 
(right-hand figure) corresponds to the initial data 

, r -2, x <o, 

MX) = { 4, *>0, 

and the corresponding solution is a classical shock followed by a nonclassical 
shock. 

Test F. We now study how the kinetic relation ur = (p (ul) is computed. 
On Figure [9] (right-hand figure) , we plot points whose horizontal coordinates 
(respectively vertical coordinates) correspond to the left-hand (resp. right- 
hand) traces around the reconstructed cell. The initial data allows us to 
cover a large range of value: 

{0, x<0.5, 
1 + 20(x + 0.45), 0.5 < x < 0.45, 
-0.75, x > -0.45. 

The left-hand figure represents the solution at different times with Ax = 
0.0002. 

We clearly observe the convergence of the numerical kinetic relation 
towards the prescribed one. This a strong test to validate the proposed 
method. 

Test G. In the course of designing the scheme proposed in the previous 
section we tried several variants. We report here one such scheme that 
is very similar to the proposed scheme, but which does not converge to 
exact nonclassical solutions. This is due to the fact that small oscillations 
are generated in the scheme which are in competition with the dissipation 
mechanisms described by the prescribed kinetic function. 

The variant is designed for the concave-convex flux f(u) = u 3 + u. 
The only difference with the scheme developed above is that it performs 
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att=0.00 - 
exact solution at t=0.01 
numerical at t=0.01 - 500 points - 



at t=0.00 - 
exact solution at t=0. 01 - 
numerical at t=0.01 - 500 points 



?ure 8: Test E - Two examples in the convex-concave case 
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the reconstruction in Cj with 



i r l_ 1 (instead of tp 



) and u 



<P (Uj_i). This is equivalent in the case of a pure nonclassical shock (Test 
B) but different in the general case. 



Figure 10 presents the solution obtained for the same initial value as in Test 
E. Oscillations are generated because the reconstruction is not constrained 
enough in this version of the scheme. 




u(0.00,x) 
u(0.05,x) - 5 000 mailles - 
u(0.05,x) - 50 000 mailles - 



0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 



Figure 10: Another version of the scheme 



6 Concluding remarks 

In this paper we have introduced a new numerical strategy for comput- 
ing nonclassical solutions to nonlinear hyperbolic conservation laws. The 
method is based on a reconstruction technique performed in each computa- 
tional cell which may exhibit a nonclassical shock. Importantly, the whole 
algorithm is conservative and propagates any admissible nonclassical dis- 
continuity exactly. The convergence of the proposed method was demon- 
strated numerically for several test-cases. This new approach brings a new 
perspective on the numerical approximation of nonclassical shocks and ki- 
netic functions. The efficiency of the method is clearly demonstrated in the 
present paper, and we refer to the follow-up paper [3] for various extensions 
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and applications. Among the questions of interest we can mention the to- 
tal variation bounds and the hyperbolic systems of conservation laws, the 
application to real materials undergoing phase transitions, as well as the 
extension to higher-order schemes. 
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